######## Attitudes to crime and punishment in England and Wales, 1964-2023
######## 004. Analysis_MediationModel

library(tidyverse)
library(gt)
library(ggraph)
library(tidygraph)
library(haven)


####### Correlation amongst attitudinal variables

Trends <- read.csv("/Users/matteo/Desktop/Prison Polling/LatentTrendsAllAdults.csv")

Trends %>% # 
  select(Punitiveness, CrimeConcern, MIP) %>%
  rename(concern = CrimeConcern, punish = Punitiveness,
         priority = MIP) %>%
  psych::corr.test(., adjust = "none", ci = TRUE) %>% 
  print(short=FALSE)


####### Mediation models

### Pooled BSA analysis

read_csv("/Users/matteo/Desktop/Prison Polling/Prison Data Internship/Survey Data/Finished/BSA_allrepeatedcrimequestions.csv") %>%
  mutate(spend1 = ifelse(spend1 %in% c('police and prisons','police & prisons'), 1, 
                         ifelse(is.na(spend1), NA_integer_, 0)),
         spend2 = ifelse(spend2 %in% c('police and prisons','police & prisons'), 1, 
                         ifelse(is.na(spend2), NA_integer_, 0)),
         govspend = spend1 + spend2,
         stifsent = ifelse(stifsent %in% c("agree", "agree   strongly","agree strongly"), 1,
                           ifelse(is.na(stifsent), NA_integer_, 0)),
         burglary = ifelse(burglary %in% c("fairly common", "very common", "fairly  common", "very    common"), 1,
                           ifelse(is.na(burglary), NA_integer_, 0)),
         rage = rage,
         rsex = case_match(rsex,
                           'Male' ~ 1,
                           'Female' ~ 0),
         Date = Year,
         ethnicgp = case_match(ethnicgp,
                               'white' ~ 0,
                               c('black', 'asian', 'other') ~ 1)) %>%
  select(Date, rsex, rage, ethnicgp, wtfactor, burglary, govspend, stifsent) -> readforanalysis

readforanalysis %>% 
  select(!c(wtfactor)) -> df1
df1[complete.cases(df1),]-> df1 # n = 5152

readforanalysis %>%
  group_by(Date) %>%
  summarise(across(where(is.numeric), \(x) weighted.mean(x, wtfactor, na.rm = TRUE))) %>%
  select(!c(wtfactor,rsex,rage,ethnicgp)) -> means

df1 %>% 
  left_join(., means, by = 'Date') %>%
  mutate(Punitiveness = stifsent.x - stifsent.y,
         MIP = govspend.x - govspend.y,
         Fear = burglary.x - burglary.y) %>%
  select(Date,rsex,rage,ethnicgp,Fear,MIP,Punitiveness) -> df1
# covers 1996 1994 1990


mediation.model <- '
  # direct effect
  Punitiveness ~ c*Fear
  # mediation
  MIP ~ a*Fear + c1*rage + c2*rsex + c3*ethnicgp
  Punitiveness ~ b*MIP + c3*rage + c4*rsex + c5*ethnicgp

  # indirect and total effects
  ab := a * b
  total := c + ab
  '

mediation.fit_BSA <- lavaan::sem(mediation.model, data=df1,
                                  se = "bootstrap", test = "bollen.stine",
                                  bootstrap = 50000)
lavaan::summary(mediation.fit_BSA, standardized=TRUE, fit.measures = TRUE)
# Appendix 11: Testing the ‘transmission mechanism’ model (BSA)


### BESIP analysis

read_csv("/Users/matteo/Desktop/Prison Polling/Prison Data Internship/Survey Data/Finished/BESIP_allrepeatedcrimequestions.csv") %>%
  mutate(gender = case_when(gender == 1 ~ 1, gender == 2 ~ 0)) %>%
  mutate(across(starts_with('mii_cat'),
                ~if_else(. == 14, 1, 
                         ifelse(is.na(.), NA_integer_, 0)), .names = "crime_{col}")) %>% # 1 if mii is crime, 0 otherwise
  mutate(across(starts_with('changeCrimeW'),
                ~case_match(., c(4,5)~1, c(1,2,9999)~0))) %>% # 1 if getting a little / a lot higher, 0 otherwise
  mutate(across(starts_with('al5'),
                ~case_match(., c(4,5)~1, c(1,2,9999)~0))) %>% # 1 if agree or agree strongly with stiffer sentences, 0 otherwise
  mutate(across(starts_with('p_ethnicityW'),
                ~case_when(. %in% c(1,2) ~ 0,
                           . %in% c(3,4,5,6,7,8,9,10,11,12,13,14,15) ~ 1))) %>%
  select(id, gender, starts_with('ageW'), starts_with('p_ethnicityW'),
         starts_with('crime_'), starts_with('changeCrimeW'),
         starts_with('al5')) %>%
  pivot_longer(cols = !c(id, gender), 
               names_to = c('variable', 'wave'),
               names_sep = "W") %>%
  mutate(wave = as.numeric(wave),
         date = case_match(wave,
                           1 ~ "2014-02-01",
                           2 ~ "2014-05-01",
                           3 ~ "2014-10-01",
                           4 ~ "2015-03-01",
                           5 ~ "2015-04-01",
                           6 ~ "2015-05-01",
                           7 ~ "2016-04-01",
                           8 ~ "2016-06-01",
                           9 ~ "2016-07-01",
                           10 ~ "2016-12-01",
                           11 ~ "2017-04-01",
                           12 ~ "2017-05-01",
                           13 ~ "2017-06-01",
                           14 ~ "2018-05-01",
                           15 ~ "2019-03-01",
                           16 ~ "2019-06-01",
                           17 ~ "2019-11-01",
                           18 ~ "2019-11-20",
                           19 ~ "2019-12-01",
                           20 ~ "2020-06-01",
                           21 ~ "2021-05-01",
                           22 ~ "2021-12-01",
                           23 ~ "2022-05-01",
                           24 ~ "2022-12-01",
                           25 ~ "2023-05-01")) %>%
  drop_na(date) %>%
  select(-wave) %>%
  pivot_wider(id_cols = c(id, gender, date), 
              names_from = variable, values_from = value) -> subset

read_csv("/Users/matteo/Desktop/Prison Polling/Prison Data Internship/Survey Data/Finished/BESIP_summaries.csv") %>%
  filter(Demographic == "All adults") %>%
  select(Date, Varname, Index) %>%
  pivot_wider(names_from = Varname, values_from = Index) %>%
  rename(crime_mii_cat = mii_crime, changeCrime = increasecrimeBSIP,
         al5 = stiffsentBESIP) %>%
  mutate(date = as.character(Date)) %>%
  select(date, crime_mii_cat, changeCrime, al5) -> means

left_join(subset, means, join_by(date)) %>%
  mutate(Punitiveness = al5.x - al5.y,
         MIP = crime_mii_cat.x - crime_mii_cat.y,
         Fear = changeCrime.x - changeCrime.y) %>%
  select(id, date,age,gender,p_ethnicity,Fear,MIP,Punitiveness) -> df2

df2[complete.cases(df2),]-> df2
# 32% of respondents answered the questions more than once
# so take average of those cases
df2 %>%
  arrange(as.Date(date)) %>%
  group_by(id) %>%
  summarise(Fear = first(Fear),
            MIP = first(MIP),
            Punitiveness = first(Punitiveness),
            age = first(age),
            gender = first(gender),
            p_ethnicity = first(p_ethnicity)) -> df2


mediation.model <- '
  # direct effect
  Punitiveness ~ c*Fear
  # mediation
  MIP ~ a*Fear + c1*age + c2*gender + c3*p_ethnicity
  Punitiveness ~ b*MIP + c3*age + c4*gender + c5*p_ethnicity

  # indirect and total effects
  ab := a * b
  total := c + ab
  '

mediation.fit_BESIP <- lavaan::sem(mediation.model, data=df2,
                                   se = "bootstrap", test = "bollen.stine",
                                   bootstrap = 50000)
lavaan::summary(mediation.fit_BESIP, standardized=TRUE, fit.measures = TRUE)
# Appendix 12: Testing the ‘transmission mechanism’ model (BES Internet Panel)



####### Prepare path model diagrams

lavaan::summary(mediation.fit_BSA, standardized=TRUE)$pe %>%
  filter(op %in% c("~","=~"), lhs != rhs) %>%
  transmute(to = lhs,
            from = rhs,
            weight = abs(std.all),
            type = case_when(op == "~" ~ 'reg', op == "=~" ~ 'load'),
            linetype = case_when(std.all > 0 ~ 'pos', std.all < 0 ~ 'neg'),
            labels = paste0(round(std.all, digits = 2),
                            ifelse(is.na(pvalue),'',
                                   ifelse(pvalue < 0.05 & pvalue >= 0.01,' *',
                                          ifelse(pvalue < 0.01 & pvalue >= 0.001, " **",
                                                 ifelse(pvalue < 0.001, " ***",'')))))) %>%
  mutate(to2 = ifelse(type == "load", from, to),
         from2 = ifelse(type == "load", to, from),
         to = to2, from = from2,
         labels = ifelse(type == 'load','',labels)) %>%
  select(-c(from2, to2)) %>%
  mutate(across(c(to,from), ~ case_match(., "Punitiveness" ~ "Punish", 
                                         'Fear' ~ 'Concern',
                                         'MIP' ~ 'Priority'))) %>%
  drop_na(from) %>%
  tbl_graph(edges = .) %>%
  activate(nodes) %>%
  create_layout(., layout = 'manual', circular = FALSE,
                x = c(0,2,1), y = c(1,1,2)) %>%
  ggraph() +
  geom_edge_link(aes(label = labels),
                 angle_calc = "along", vjust = -0.5,
                 arrow = arrow(length = unit(4, "mm")),
                 end_cap = circle(10, 'mm')) +
  geom_node_point(shape = 22, fill = "black", size = 20, color = "black") +
  geom_node_text(aes(label = name), color = 'white', size = 3) +
  scale_y_continuous(limits = c(0,2.5)) +
  scale_x_continuous(limits = c(-1,3)) +
  theme_void() +
  labs(title = "A. Individual level analysis, BSA (1990-1996, n = 5143)",
       caption = "Notes: CFI = 0.974, TLI = 0.765, RMSEA = 0.029.\nControls for gender, age, race not shown.") +
  theme(legend.position = "none",  plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, vjust = 20)) -> p1

lavaan::summary(mediation.fit_BESIP, standardized=TRUE)$pe %>%
  filter(op %in% c("~","=~"), lhs != rhs) %>%
  transmute(to = lhs,
            from = rhs,
            weight = abs(std.all),
            type = case_when(op == "~" ~ 'reg', op == "=~" ~ 'load'),
            linetype = case_when(std.all > 0 ~ 'pos', std.all < 0 ~ 'neg'),
            labels = paste0(round(std.all, digits = 2),
                            ifelse(is.na(pvalue),'',
                                   ifelse(pvalue < 0.05 & pvalue >= 0.01,' *',
                                          ifelse(pvalue < 0.01 & pvalue >= 0.001, " **",
                                                 ifelse(pvalue < 0.001, " ***",'')))))) %>%
  mutate(to2 = ifelse(type == "load", from, to),
         from2 = ifelse(type == "load", to, from),
         to = to2, from = from2,
         labels = ifelse(type == 'load','',labels)) %>%
  select(-c(from2, to2)) %>%
  mutate(across(c(to,from), ~ case_match(., "Punitiveness" ~ "Punish", 
                                         'Fear' ~ 'Concern',
                                         'MIP' ~ 'Priority'))) %>%
  drop_na(from) %>%
  tbl_graph(edges = .) %>%
  activate(nodes) %>%
  create_layout(., layout = 'manual', circular = FALSE,
                x = c(0,2,1), y = c(1,1,2)) %>%
  ggraph() +
  geom_edge_link(aes(label = labels),
                 angle_calc = "along", vjust = -0.5,
                 arrow = arrow(length = unit(4, "mm")),
                 end_cap = circle(10, 'mm')) +
  geom_node_point(shape = 22, fill = "black", size = 20, color = "black") +
  geom_node_text(aes(label = name), color = 'white', size = 3) +
  scale_y_continuous(limits = c(0,2.5)) +
  scale_x_continuous(limits = c(-1,3)) +
  theme_void() +
  labs(title = "B. Individual level analysis, BESIP (2018-2023, n = 40290)",
       caption = "Notes: CFI = 0.996, TLI = 0.967, RMSEA = 0.019.\nControls for gender, age, race not shown.") +
  theme(legend.position = "none",  plot.title = element_text(hjust = 0.5),
        plot.caption = element_text(hjust = 0.5, vjust = 20)) -> p2

rm(df1,df2,lmready,means,readforanalysis,subset,Trends,mediation.model)


### Figure 7: Path diagrams for two tests of the ‘transmission mechanism’ hypothesis

cowplot::plot_grid(p1,p2,ncol=1)


save.image(file = "/Users/matteo/Desktop/Prison Polling/MediationModels.RData")
